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We calculate superconducting transition temperatures (T c ) in sulfur hydrides H 2 S and H 3 S from 
first principles using the density functional theory for superconductors. At pressures of < 150 GPa, 
the high values of T c (>130 K) observed in the recent experiment [A. P. Drozdov, M. I. Eremets, 
and I. A. Troyan, arXiv: 1412.0460 are accurately reproduced by assuming that H 2 S decomposes 
into i? 3 m-H 3 S and S. For the higher pressures, the calculated T c s for 7 m 3 m-H 3 S are systematically 
higher than those for J? 3 m-H 3 S and the experimentally observed maximum value (190 K), which 
suggests the possibility of another higher- T c phase. We also quantify the isotope effect from first 
principles and demonstrate that the isotope effect coefficient can be larger than the conventional 
value (0.5) when multiple structural phases energetically compete. 

PACS numbers: 74.62.Fj, 74.20.-z, 74.25.Kc, 74.70.-b 


I. INTRODUCTION 

Investigating compounds containing light elements has 
been a simple and powerful guiding principle for discov¬ 
ery of high-temperature superconductors. According to 
the BCS theoryfi the superconducting transition temper¬ 
ature ( T c ) is scaled by the phonon frequency and there¬ 
fore light atoms are advantageous for achieving high T c . 
Despite its simplicity, this principle has been surpris¬ 
ingly successful as represented by the discoveries of su¬ 
perconductivity in doped fullerene solids^ magnesium di- 
boride^ lithium under pressure^ and boron-doped dia¬ 
mond^’' Along this principle, possible superconductivity 
in compressed hydrogen and hydrogen compounds has 
been explored as an extreme case£~— 

Recently, it has been discovered that H 2 S exhibits su¬ 
perconductivity under high pressures at 190K (Ref. [36). 
Being the new record of the superconducting transition 
temperature (T c ), this report has immediately aroused 
intense debate^^ Several facts imply that this super¬ 
conducting phase is induced by the conventional mecha¬ 
nism due to the vibrations of hydrogen atoms: The ob¬ 
served T c is subject to the hydrogen isotope effect^; in 
prior to the experimental discovery, there was an ab ini¬ 
tio calculation which predicted strong electron-phonon 
coupling^!; the electronic bandwidth is so large that the 
Migdal approximation seems valid4 2 However, some puz¬ 
zling results have also been exposed. First, the crystal 
structure realized in the experimental situation has not 
been specified. If we estimate T c of H 2 S using the conven¬ 
tional McMillan formul a 34 : 43 with the empirical Coulomb 
parameter p*=0.13, the calculated value is too low com¬ 


pared with the experimentally observed value. It has 
also been proposed that H 3 S phase instead emerges under 
high pressures , 39 ' 41 where the electron-phonon coupling 
is thought to be stronger than in Second, anoma¬ 

lously large hydrogen isotope effect coefficient a~ 1.0 has 
been observed. Although it has been hypothesized that 
the unharmonic effect on the lattice dynamics has some 
role^i or that different structures emerge in H 2 S and D 2 S 
(sulfur di-deuteride)j 2 £ this anomaly remains an open 
question. 

To further investigate the above points, we need to ad¬ 
dress not only the electron-phonon interaction but also 
the electron-electron Coulomb interaction in the H^S sys¬ 
tems. Accurate evaluation of the impact of the pair¬ 
breaking Coulomb repulsion is vital because this governs 
the absolute value of T c , as well as a: 44- — In addition, 
experimentally realized pressure range is rather out of 
that in the previous ab initio studies and therefore more 
thorough investigations of the pressure dependence of the 
superconducting properties are desired. 

In this Article, we present an ab initio study on the 
superconductivity in solid H 2 S and H 3 S covering the 
experimental pressure range. In the standard Migdal- 
Eliashberg theory ; 42 : 47 the effect of the electron-electron 
Coulomb interaction is practically treated with an em¬ 
pirical parameter /x*. To incorporate this effect non- 
empirically, we utilize the density functional theory for 
superconductors (SCDF T 48 : 49 ). With this theory, we 
can calculate T c and a without any empirical parameter, 
which can be directly compared with the experimental 
data. 
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II. METHOD 


To calculate T c from first principles, we employed the 
SCDFT gap equation given by 


Here, n and k denote the band index and crystal mo¬ 
mentum, respectively, A„k is the gap function, and /3 
is the inve rse tempera ture. The energy E n k is defined 
as ^«k = V£nk + A nk and £nk is the one-electron en¬ 
ergy with respect to the Fermi level calculated with the 
normal Kohn-Sham equation. The functions Z and JC 
are called exchange-correlation kernels, which describe 
the effects of the interactions. The nondiagonal kernel 
K. consists of two parts /C=/C ph +/C el representing the 
electron-phonon and electron-electron interactions, re¬ 
spectively. The diagonal kernel Z=Z ph represents the 
mass renormalization of the normal-state band structure 
due to the phonon exchange. Using these kernels ; 48 ’ 49 
the conventional strong-coupling superconductivity can 
be treated with the level of the Migdal-Eliashberg the- 
oryi 42 ’ 47 In particular, the electronic nondiagonal kernel 
/C el describes the screened electron-electron Coulomb in¬ 
teraction, where the dynamical screening effects are in¬ 
corporated within the random-phase approximatio m 50 ’ 51 
We can therefore evaluate effects of the static Coulomb 
repulsion suppressing the pairing, as well as the plasmon 
superconducting mechanism— 

We calculated the electronic states, phonon frequen¬ 
cies, electron-phonon and electron-electron interactions 
and T c for H 2 S and H 3 S at various pressures. Our cal¬ 
culations were performed with the generalized-gradient 
approximation using the exchange-correlation poten¬ 
tial with the Perdew-Burke-Ernzerhof parametrizatiom 53 
We used ab initio plane-wave pseudopotential calcula¬ 
tion codes QUANTUM ESPRESSO^ for the electroic 
structure, dynamical matrix and electron-phonon cou¬ 
pling. The input crystal structures at respective pres¬ 
sures were the optimum ones predicted in the previous 
ab initio calculations, which are summarized in Table [U 
For all the conditions, we optimized the atomic configura¬ 
tions and cell parameters with respect to enthalpy under 
fixed pressures. Phonon frequencies and electron-phonon 
interactions were calculated based on the density func¬ 
tional perturbation theory— The electron dielectric func¬ 


TABLE I: Pressure settings and the corresponding input 
structures for the calculations. We observed that it is difficult 
to achieve the numerical convergence in the phonon calcula¬ 
tions for the calculations for H 3 S at 190 GPa since it is near 
the second-order structural transition point — 



FIG. 1: Calculated superconducting transition temperatures 
for H 2 S (solid square) and H 3 S (solid circle). Experimen¬ 
tally observed values for H 2 S [Fig. 2 (a) (open square) and 
Fig. 2(b) (open circle) of Ref. |3f| are also plotted together, 
where different runs are represented by the same symbols. 
Open pentagon and diamond denote the ab initio predictions 
for H 2 £— and H 3 S,— respectively. The small solid circle for 
the 7 m 3 m-H 3 S phase indicates the calculated result without 
the contribution of the plasmon mechanism. 


tions were calculated within the random-phase approx¬ 
imation, where the frequency dependence was retained. 
K. ph and Z ph were calculated with the nk-averaged ap¬ 
proximate formula [Eq. (23) in Ref. |49| and Eq. (40) in 
Ref. !H, respectively], whereas /C el was calculated includ¬ 
ing the plasmon-induced dynamical screening effec t 50 ’ 51 . 
The SCDFT gap equation was solved with the random 
sampling scheme given in Ref. l57l with which the sam¬ 
pling error was approximately a few %. Further details 
are summarized in Appendix lAl 

We took particular care for calculating the Eliashberg 
function 

nn'k 

employed for /C ph and Z ph . N( 0), 9 ”£ +q ’”' k and u„ q 
denote the density of states at the Fermi energy, the 
electron-phonon matrix element and the phonon fre¬ 
quency, respectively. Since we have found that a 2 F(u) 
sensitively depends on the smearing scheme and k- and q- 
point density for the integration, we employed a recently 
developed tetrahedron method with optimized linear in¬ 
terpolation^ 

We included the plasmon-induced frequency depen¬ 
dence of the screened Coulomb interaction in K. el with 
the following formula [Eq. (2) of Ref. [50j ] 


P [GPa] 

130|150 

170 | 190 | 210 

230 | 250 

H 2 S 

Pl^ 

Cmca f 


h 3 s 

R3m^ 


Im3m — 


ir-el,dvn ,. 1 1 

p, ' j — ii m _ 

nk,n'k A„ k ->o tanh[(/3/2).E„k] tanh[(/3/2)£ , ri 
x ^2 E rl k(iWl)F„k(iW 2 )H ; nkrt'k'[i(wi^U 2 )], 


(3) 


UJ1UJ2 
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where W^kn'k' [i(wi—i^ 2 )] is the screened Coulomb inter¬ 
action and F n k(iLo)= - -^e — denote the elec¬ 

tronic anomalous Green’s function, m the previous cal¬ 
culations we carried out the Matsubara summations 
analytically by approximating W n kn'k'[i( w i — w 2 )] with 
model functions. In the present study, the summation 
for ui\ was done analytically with variable transformation 
u>± — w 2 = v, whereas the summation for v was evaluated 
numerically with ~ y / di/ without any modeling of 
bCnkn'k'[i(wi—w 2 )], where T is the temperature— 


III. RESULTS AND DISCUSSION 

Below, we show the calculated values of T c and key 
factors for the phonon theory: A, wi n , /a* and the isotope 
effect coefficient a. The specific values are summarized 
in Appendix [Bl 

In Fig. [[I we show the calculated T c with the previ¬ 
ously published experimental and first-principles numeri¬ 
cal data .— ~- 6 Drozdov and coworkers^ reported two data 
groups obtained with different experimental conditions, 
which are indicated by open square and circle, respec¬ 
tively; in this work, we name these groups data 1 and 
data 2, respectively. The calculated T c s for H 2 S (solid 
square) and H 3 S (solid circle) were —50 K and >130 K, 
respectively. For both H 2 S and H 3 S, the calculated T c s 
show domelike dependence on the pressure. The max¬ 
imum T c s are achieved near the theoretically proposed 
structural transition points i 34 ’ 35 Our calculated values 
are as a whole in good agreement with the previous esti¬ 
mates with the McMillan-Alien-Dynes formula (Refs. l34l 
and [351. Notably, for H 3 S, we obtained 267 K at max¬ 
imum, which is larger by —60 K from the previous es¬ 
timate^ at 200 GPa. This difference is discussed more 
specifically later. In the low-pressure regime, the cal¬ 
culated T c for H 2 S (H 3 S) agrees well with data 1 (data 
2). In the high-pressure regime, on the other hand, the 
calculated values are too high or too low compared with 
the experimental ones. Furthermore, the rapidly increas¬ 
ing feature of data 1 (> 170 GPa) was not reproduced. 
We also revisit this point later. Regarding the plasmon 
effect ; 50 ! 51 the enhancement of T c was estimated to be 
15-20% (-10%) for H 2 S (H 3 S) (e.g., see small solid cir¬ 
cle). 

To understand the pressure dependence of the calcu¬ 
lated T c in terms of the McMillan-Alien-Dynes formula^ 
we show the calculated values of A = 2 f duj a F ^ and 

wi n = exp[j / duj a F fcd lrm; ] i n Figs. [2] (a) and (b), re¬ 
spectively. We see that the pressure dependence of the 
present ab initio T c s for H 2 S and H 3 S are similar to that 
of A, which indicates that the T c s of the present sys¬ 
tems are governed by A. With this plot for A, we see 
that our tetrahedron method^ and the previously em¬ 
ployed Gaussian smearing schem o 35 ! 60 give different val¬ 
ues for A, which results in the large difference in T c . 
In fact, by calculating A with the first-order Hermite- 
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FIG. 2: Key factors in the conventional theory for the phonon 
mechanism calculated from first principles: (a) A, (b) wi n , (c) 
n* and (d) a. Solid square (circle) denote the values for H 2 S 
(H 3 S). Open pentagon and diamond represent the preceding 
ab initio calculations for FBS^ and H 3 S,— respectively. 
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FIG. 3: Numerical convergence of A with different schemes for 
the phonon and a 2 F(u) calculations: Optimized tetrahedron 
and the lst-order Hermite-Gaussian smearing with width of 
0.030 Ry. kph and k ep represent the k-point grids employed 
for the phonon dynamical matrix and Eliashberg function, re¬ 
spectively. The q-point summation for “Smearing” was done 
with a q-point grid without offset. 
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Gaussian approximate function [£(£) ~ [3/2 — 

(^/VF) 2 ]exp(—(^/VF) 2 ) with IV=0.030 Ry, Ref. IH3], we 
obtained A=2.23 and 1.99 for P =200 GPa and 210 GPa, 
respectively, which is consistent with the previous value 
(A=2.19 for P=200 GPa—). Since the bandwidth of 
the electronic states is extremely large and complex- 
shaped electron/hole pockets emerge in this system^ 
the present tetrahedron-interpolation-based method is 
expected to be more numerically accurate. We have 
confirmed the numerical convergence of A as depicted 
in Fig. [3] Wi n monotonically increases as the pressure 
is increased, which represents the hardening of phonons 
by compression. This hardening is responsible for the 
marked difference in T c s for P 3771 -H 3 S and ImSm-UsS. 
For higher pressure regime, however, the hardening is 
dominated over by the decrease of A and therefore T c 
decreases— 

We determined optimum values for /i* so that the T c s 
calculated with the SCDFT gap equation can be repro¬ 
duced with the extended McMillan formula— For H 2 S, 
the optimum values were 0.15-0.17 for all the pressures. 
For the pressure range 170-210 GPa, we observed de¬ 
crease of the optimum values for H 3 S. Probably this is 
originating from a fact that T c calculated by the present 
SCDFT sometimes deviates slightly from that by the 
Eliashberg equation— Detailed investigations on this 
point are left for future studies. 

Using the calculated T c s for H X S and D X S, we also 
calculated the isotope effect coefficient a = — [lnT c DxS — 
lnT c HxS ]/[lnMD — lnMn], where T c HxS (T c DxS ) is the tran¬ 
sition temperature of the hydride (deuteride) compound 
and Mh (Mq) is the atomic mass of hydrogen (deu¬ 
terium), respectively. The values ranges between 0.23 
and 0.31 (0.38 and 0.42) for H 2 S (H 3 S). These values are 
smaller than the BCS value (a^0.5), which indicate the 
correction due to the retardation effect. 

Here we compare our calculated and experimentally 
observed values of T c . First, the experimentally observed 
T c s in the low-pressure regime were quantitatively re¬ 
produced by assuming the emergence of single structural 
phases of P1-H 2 S and RSm-^S for data 1 and 2, respec¬ 
tively. This strongly suggests that these two phases are 
dominant in the experimental situations for P<150 GPa. 
It is even conceivable that the high-pressure values of 
data 2 corresponds to P3m-H 3 S. The agreement of the 
calculated and experimentally observed T c s for higher 
pressures were, on the other hand, not as perfect as 
those for the previously studied conventional supercon¬ 
duct or s42i5Si§Irij£ Note that we assumed that the sample 
is homogeneous and does not decompose into H X S and 
S for all the pressure range, though it has not been con¬ 
firmed experimentally. Our calculated T c for Im3m-H 3 S 
suggests that maximum T c can be increased to, possibly, 
a higher value in the pure 7 ro 3 m-H 3 S phase. 

Very recently, there has been an independent report 
on an ab initio T c calculation for Im3m-H 3 S using the 
SCDFT— with a condition different from ours— They 
concluded that the experimentally observed high T c can 


be explained with 7 m 3 ro-H 3 S, whereas we rather propose 
a relevance of R3m-H 3 S in the experimental situation. 

Finally, we move on to a. The calculated values 
were far smaller than the experimentally observed a~ 1 . 0 . 
Based on a hypothesis of inhomogeneity, let us give a 
possible explanation of the experimental large a within 
the present framework. As suggested by Hirsch and 
MarsigliOf?- when inhomogeneity of the system is sub¬ 
stantial, the experimentally observed T c should some¬ 
how deviate. For example, suppose we estimate a with 
a = — [lnTjP 2S — lnT c H3S ]/[lnMD — lnMn]; we then get 
a> 2.0 for the whole pressure range. Such a situation is 
possible because the enthalpy difference between H 2 S and 
§H 3 S+±S is of order of the phonon frequency— Substi¬ 
tution of D for H substantially modulate the contribution 
of the zero-point oscillation to the total enthalpy and 
it should change the relative stability of the competing 
phases. 

We thus suggest that the H 3 S phases have a key role 
for understanding the reported experimental results^ 
and realizing higher T c . To validate/invalidate this, 
measurements with different chemical composition (e.g., 
H:S=3:1) and compression at higher temperatures might 
be helpful. When measuring the isotope effect, the dif¬ 
ference in the structural relaxation speed of hydrides and 
deuterides should also be taken into account. 


IV. SUMMARY 


In this study, we have performed a present state-of- 
the-art ab initio calculation for the superconductivity in 
H 2 S and H 3 S assuming the conventional phonon mecha¬ 
nism, where the effect of the electron-electron Coulomb 
repulsion was non-empirically treated. For the pres¬ 
sures < 150 GPa, the calculated T c s for P1-H 2 S and 
-R 3 TO-H 3 S agree well with the experimental T c s observed 
with different compressing and cooling conditions, re¬ 
spectively. This strengthens the scenario that H 3 S is 
superconducting when the high T c is observed— ^ For 
the high-pressure phase of Jm 3 m-H 3 S, we have predicted 
T c higher than the experimentally observed maximum 
of 190 K and the values calculated for R3m-H 3 S, which 
amounts to 267 K. This suggests that higher T c can be 
achieved by isolating the single Im3m-H 3 S phase. Al¬ 
though we have ignored several possible effects in the 
present systems (e.g., zero-point oscillation of hydrogen 
atoms, anharmonic phonons, etc.), the present result can 
be a key step for further theoretical and experimental 
investigations on the superconducting sulfur hydrides. 
Examinations of anharmonic lattice-dynamical effects, 
which has been neglected with the present methodology, 
are under way. 
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TABLE II: Detailed settings for the calculations. Subscript “1” for q points denotes the mesh with displacement by half a grid 
step. 




P 1-H 2 S 

Cmca- H 2 S 

R3m-H 3 S 

Im3m H 3 S 

charge density 

k 

(12 12 8 ) 

(12 12 4) 

(16 16 16) 

(16 16 16) 

interpol. 

1 SI order Hermine Gaussian^ with width=0.030Ry 

dynamical matrix 

k 

(12 12 8 ) 

(12 12 4) 

(16 16 16) 

(16 16 16) 

q 

(6 6 4)i 

(6 6 2 )r 

(8 8 8 )i 

(8 8 8)1 

interpol. 

Optimized tetrahedron^ 

electron-phonon 

k* 

(12 12 8) 

(24 24 8) 

(32 32 32) 

(32 32 32) 

interpol. 

Optimized tetrahedron^ 

dielectric function 

k for bands crossing Ep 

(18 18 12) 

(18 18 6) 

(18 18 18) 

(18 18 18) 

k for other bands 

(6 6 4) 

(6 6 2) 

(6 6 6) 

(6 6 6) 

q 

(6 6 4) 

(6 6 2) 

(6 6 6) 

(6 6 6) 

unoccupied band num. 

-60 

-100 

—30 

-30 

interpol. 

Tetrahedron with the Rath-Freeman treatment 

SCDFT gap function 

unoccupied band num. 

25 

45 

19 

19 

k for 1C' 

(6 6 4) 

(6 6 2) 

(6 6 6) 

(6 6 6) 

N s for bands crossing Ep 

4500 

3000 

6000 

6000 

N s for other bands 

150 

100 

200 

200 

Sampling error in T c 

-9% 

-6% 

-5% 

-5% 


t Electron energy eigenvalues and eigenfunctions were calculated on these auxiliary grid points. 
^ Electron energy eigenvalues were calculated on these auxiliary grid points. 


Note added 

After the submission of the present work, there has 
been a publication demonstrating that the anharmonic 
effect reduces T c by about 20 % in /m3m-H 3 S^ Nev¬ 
ertheless, the present indications of the relevance of the 
R3m phase and possible higher T c are still valid since the 
increase of T c by the plasmon mechanism will compen¬ 
sate the anharmonic effect. 
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Appendix A: Computational detail 

For the electronic and lattice-dynamical calculations, 
we used the pseudopotentials for S and H atoms imple¬ 
mented with the Troullier-Martin scheme^- which are 
the same as those used in Ref. [HI- The plane-wave en¬ 
ergy cutoff was set to 80 Ry, whereas the auxiliary cutoff 
for the dielectric function was 12.8 Ry. Conditions for 
the calculations of the charge density, dynamical ma¬ 
trix, electron-phonon coupling, dielectric function and 
gap function are detailed in Table HT1 


Appendix B: Numerical data of the calculated 
values of T c , A, win, p* and a 

Tables HTT1 1 VIII lists the calculated values for T c , A, win, 
/r and a. 

TABLE III: Superconducting transition temperature T c [K]. 


P [GPa] 

130 

150 

170 

190 

210 

230 

250 

H 2 S 

29.4 

47.1 

66.9 

56.3 

55.4 

45.4 

41.8 

d 2 s 

25.0 

38.2 

55.7 

46.7 

45.0 

37.9 

35.7 

h 3 s 

134 

155 

214 


267 

236 

211 

d 3 s 

103 

119 

163 


206 

180 

164 


TABLE IV: Electron-phonon coupling coefficient A. 


P [GPa] 

130 

150 

170 

190 

210 

230 

250 

H 2 S 

0.801 

1.001 

1.196 

1.026 

0.945 

0.882 

0.855 

h 3 s 

1.843 

2.067 

2.599 


2.582 

2.263 

1.970 


TABLE V: Logarithmic moment of the Eliashberg function 
Win [K], 


P [GPa] 

130 

150 

170 

190 

210 

230 

250 

H 2 S 

913 

914 

968 

1044 

1097 

1121 

1124 

h 3 s 

1037 

1056 

1058 


1336 

1447 

1521 
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TABLE VI: Renormalized electron-electron Coulomb param¬ 
eter estimated from the T c calculated with the SCDFT 
gap equation. 


P [GPa] 

130 

150 

170 

190 

210 

230 

250 

H 2 S 

0.155 

0.165 

0.174 

0.166 

0.152 

0.158 

0.159 

h 3 s 

0.168 

0.165 

0.125 


0.118 

0.153 

0.164 


TABLE VII: Isotope-effect coefficient a. 


P [GPa] 

130 

150 

170 

190 

210 

230 

250 

H 2 S 

0.24 

0.31 

0.27 

0.27 

0.30 

0.26 

0.23 

h 3 s 

0.38 

0.39 

0.40 


0.38 

0.40 

0.37 
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